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Abstract 

We propose a very simple but 'realistic' model of amphiphilic bilayers, simple enough to be able 
to include a large number of molecules in the sample, but nevertheless detailed enough to include 
molecular charge distributions, flexible amphiphilic molecules and a reliable model of water. All 
these parameters are essential in a nanoscopic scale study of intermolecular and long range electro- 
static interactions. We also propose a novel, simple and more accurate macroscopic electrostatic 
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'^ ' field for model bilayers. This model goes beyond the total dipole moment of the sample, which on 

O , a time average is zero for this type of symmetrical samples, i. e., it includes higher order moments 

of this macroscopic electric field. We show that by representing it with a superposition of gaus- 
K*- ' sians it can be analytically integrated, and therefore its calculation is easily implemented in a MD 

^«0 . simulation (even in simulations of non-symmetrical bi- or multi-layers). In this paper we test our 
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model by molecular dynamics simulations of Newton black films. 
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I. INTRODUCTION 

The structure and dynamics of amphiphilic bilayers plays a key role in numerous 
problems of interest in physical-chemistry [1], biology [2], pharmacology [3] and biotech- 
nology [4, 5]. Amphiphilic molecules consist of a non-polar hydrophobic flexible chain of 
the {CH2)n type, the 'tail', plus a hydrophilic section, a strongly polar 'head' group. In 
aqueous solutions the polar 'head' strongly interacts with water and shields the hydrophobic 
tails from the water, therefore these molecules tend to nucleate in miscelles or bilayers, 
depending on concentration [1]. For example, a simple model of a biological membrane 
consists of a bilayer of amphiphilic molecules, with their polar heads oriented to the 
outside of the bilayer and strongly interacting with the surrounding water. The opposite 
model of bilayer, with the water in the middle and head groups pointing to the inside, 
also exists in nature and are the soap bubbles fllms, or Newton black (NB) fllms, as 
they are called in their thinnest states. Most studies have been performed on biological 
membrane models and huge advances in the intermolecular potential models (at atomistic 
or coarse grained level) needed for their numerical simulations as well as in the behavior 
of proteins, anesthetics molecules , pores, etc., inserted in these bilayers have been achieved . 

Here we will limit our study to the case of soap bubbles fllms. While on one hand 
the molecular interactions and macroscopic electrostatic held exhibit the same physical 
problems to solve as that of the biological membranes, on the other hand the amount of 
water is limited by the width of the NB fllm and therefore the numerical simulations are 
considerable less expensive. 

Furthermore, the study of foam fllms is not only interesting per se, but also has many 
technological applications, in particular their use in biothechnological applications is a very 
active held. For example, a protein can be unfolded when inserted in these fllms, allowing 
a small angle X ray measurement of its structure. In refs. [6, 7] the structure of the fllm 
is studied as a function of the concentration of inserted proteins; these measurements give 
information on the protein-lipids interactions, in much simpler experimental arrays than 
biological membranes. At the mesoscopic scale, the diffusion of bacteria in this quasi-two- 
dimensional liquid environment has been measured [8]. 



From the point of view of applied physics there is also incresing interest in the grow of 
inorganic nano-patterned films with the help of organic bilayers [9]. 

Furthermore, recent studies have revived the interest in the basic physics involved in 
these films. The self-assembly of amphiphilics can be experimentally measured by studying 
the properties of soap films, with the advantage of a fine control in the range of low contents 
and small concentration changes of amphiphilics [10]. Also confined water is known to 
have a very different behaviour of bulk water, its behavior at criogenic temperatures and 
its glass transition temperature were measured by anelastic spectroscopy [11] in a stack 
of biological membranes, but not yet in foam films. In ref. [12] the stability of soap films 
under different applied capillary pressures has been meassured, and in ref. [13] a mean field 
theory is developed to describe the role of electrostatic fluctuations in their stability. 

Several MD simulations of NB fllms have been performed since some time ago. In 
Ref. [14] a series of all-atom MD simulations of NB fllms of sodium-dodecyl-sulfate SDS 
[Na~'~CH3(CH2)ii(OS03)"] amphiphilic molecules and different amounts of water was 
presented, and an insight of the experimental electron density proflle given by X-ray 
measurements was obtained [15]. In Ref. [16] these simulations were extended to larger 
samples and longer times (about nanosec), a fact that is important to accurately measure 
the diffusion constants of the lipids in bilayers, even when similar structural and internal 
dynamical properties can be obtained at shorter times [17]. The effect of several surfac- 
tants on the disjoining pressure [18] and the anomalous behavior of water [19] were also 
investigated via MD simulations. The last reference clearly pointed out that the study of 
electrostatic forces between charged interfaces in aqueous media is still an open and very 
important fleld of theoretical research. 

Here we are interested in the search of a very simple but 'realistic' model of the Newton 
black fllms, simple enough to be able to include a large number of molecules in the 
sample, but nevertheless detailed enough to include molecular charge distributions, flexible 
amphiphilic molecules and a reliable model of water. All these parameters are essential in 
a nanoscopic scale study of intermolecular and long range electrostatic interactions. Such 
simple amphiphilic bilayer model will be useful to obtain reliable information on the effect 
of the "external parameters" (like surface tension, external pressure and temperature and 



external fields) on physical properties of the bilayer, as well as to address problems like the 
diff'usion and/or nucleation of guest molecules of technological, pharmaceutical or ambiental 
relevance within amphiphilic bilayers. 

A problem to solve in the simulation of these highly charged bilayers is the method 
to use in the calculation of electrostatic interactions, due to their long range and quasi- 
bidimensional 2D nature. The Ewald method is assumed to be the accurate one for 
electrostatic interactions [20] in three dimensional 3D samples. Refs. [21-23] review the 
3D Ewald sums as well as several possible implementations, in those papers 3D periodic 
boundary conditions are used to minimize surface effects and all interacting molecules 
are explicity included. However, for single monolayers or bilayers the periodic boundary 
conditions should be applied in two directions, in the plane of the mono- or bilayers ( x y 
plane), but not along the perpendicular to the plane {z axis). 2D Ewald sums have been 
developed for these cases, although they are computationally very lenghly by comparison 
with the 3D sums [24, 25]. Recently Brodka and Grzybowski [26] showed analytically 
in which way the 2D Ewald sums may be accurately approximated by a 3D calculation. 
They show that, in order to obtain a reliable calculation, not only a large empty space 
must be introduced in the simulation box (along z), but the macroscopic electric field term 
depending on the total dipole moment M of the sample (|M| /3) should be replaced by 
a term containing only the component z perpendicular to the 2D system (M|). In ref. 
[27] this result was checked with a numerical simulation. It has to be taken into account 
that, in a MD simulation, the contribution of this macroscopic electric field to the molec- 
ular interactions fiuctuates in time and is far from negligible [27], particularly for monolayers. 

The core of the paper is in Section IV where we propose a novel, simple and more 
accurate macroscopic field for model bilayers. That is, a model that goes beyond the total 
dipole moment of the sample, which on time average is zero for this type of symmetrical 
bilayers. Following this approach, we extend it by including higher order moments of this 
macroscopic electric field. We show that by representing it with a superposition of gaussians 
it can be analytically integrated, and therefore its calculation easily implemented in a MD 
simulation (even in simulations of non-symmetrical bilayers). 



In this paper we study a simple model of Newton black films via MD simulations and 
using the TIP5P intermolecular potential for water and a fiexible charged chain to simulate 
the amphiphilic molecules. We show that our method of calcutation of the macroscopic 
electrostatic potential reproduces the spatial and temporal charge inhomogeneities of the 
quasibidimensional sample and is extremely simple to implement in numerical simulations. 
We also show that this method can be directly applied to any other type of bi- or multi-layers. 

II. THE AMPHIPHILIC AND WATER MOLECULE MODELS 

Our simple model of a charged amphiphilic consists in a semifiexible single chain of 14 
atoms, the bond lengths are considered constant, but bending and torsional potentials are 
included. The first two atoms of the chain mimic a charged and polar head (atom 1 with qi 
=- 2 e, atom 2 with q2 = 1 e), sites 3 to 14 form the hydrophobic tail with uncharged atoms: 
sites 3 to 13 are uncharged united atom sites CH2 and site 14 is the uncharged united atom 
site CH3. This molecular model correspond to an oversimplification of sodium dodecyl 
sulfate SDS {CH2,{CH2)iiOSO^ Na'^) in solution, so we are including, in our simulations, 
a Na'^ ion per chain. The LJ parameters of these interaction sites are those of ref. [14], 
except for the sites 1 and 2 that form the amphiphilic polar head: o"i = 4.0A, o"2 =4.0 A, 
a^a =1.897 A, ei= 2.20 kJ/mol, 62= 1.80 kJ/mol and eNa= 6.721 kJ/mol. The masses of 
the sites are the corresponding atomic masses, except that mi = m2 = 48au., in order to 
mimic the 'real' amphiphilic head. The LJ parameters of the united atom sites are taken 
from calcutations on n-alkanes [28]: (Jch2 =3.850 A, acu'i =3.850 A, ecH2= 0.664 kJ/mol, 
£cH3= 0.997 kJ/mol. The LJ parameters for the Na~^ ion are taken from simulations of 
SDS in aqueous solution [29] and NB films [14]. 

The intramolecular potential includes harmonic wells for the bending angles /? and the 
usual triple well for the torsional angles r, the constants are those commonly used for the 
united atom site Ci/2[20]. These potentials are needed to maintain the amphiphilic stiffness 
and avoid molecular collapse. The bending potential is 

V{f3) = kccc{l3 - Pof , 
with Po = 109.5 deg. and kccc = 520 kJ/ra(f. The torsional potential is of the form 



V{t) = flo + aicosi^r) + a2Cos'^{T) + ascos^i^r) + a4Cos'^{T) + a^cos^^r) , 

the constants are oq = 9.2789, ai = 12.1557, aa = -13.1202, og = -3.0597, a^ = 26.2403 
and as = —31.4950 A; J/mo/; this potential has a main minimum at r = Odeg. and two 
secondary minima at r = ±120 deg. 

The selected molecular model for water is the classical rigid molecular model TIP5P 
[30, 31]. It consists of one LJ site {e = 0.67 kJ/mol, a = 3.12 A) localized at the O and 
4 charges. Two charges qpj = 0.241e are localized at the H atoms and two qip =-qH ^^ 
the lone pairs. As the lone pair interaction sites are not coincident with atomic sites, the 
algorithm employed to translate the forces from massless to massive sites is that of ref. [32]. 
The final version of the MD program is similar to that used in refs. [33-35]. 

This rigid, nonpolarizable, TIP5P molecular model was selected by us because it is simple 
and it also gives good results for the calculated energies, diffusion coeficient and density p 
as a function of temperature, including the anomaly of the density near 4C and latm [30], 
the X-ray scattering of liquid water [36], etc. The only exception is the 0-0 pair correlation 
function g2{r), for which the first neighbor is located at a slightly shorter distance than the 
experimental one [30]. Very recently a six sites rigid model for water [37], with an additional 
site at the molecular center of mass, has been presented. By comparison with the TIP5P 
model, the new model has improved the fit of the melting point and the disordered structure 
of ice at the melting point. Nevertheless, this six sites model still shows the same problem 
of TIP5P to reproduce the experimental 5'oo(^) P^^'^ correlation function and this is the 
reason to perform our calculations with TIP5P. 

III. THE MODEL BILAYER OF AMPHIPHILIC MOLECULES: 

Fig. 1 shows the initial configuration of one of our simple model of a NB film (i. e. soap 
bubbles films) with the water within the bilayer and, therefore, the head groups oriented 
to the inside. The sample of the figure includes 226 amphiphilics and 365 water molecules, 
and, as in model (a), the bilayer is perpendicular to the 'z MD box axis, with a box size of 
a = 6 = 45.A, c=1000 A. 
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Figure 1: Initial configurations of a model Newton black film. 

The periodic boundary conditions are applied only in the x and y direction and a large 
unit cell c parameter (that is, a large empty volume) is taken, in order to approximate 
the required 2D Ewald sums by the usual 3D sums plus our macroscopic field term 
corresponding to a quasi-2D system, discused in the following section. 

When performing these simulations, we have to take into account, in particular, the 
following facts: 

a) The site-site LJ interactions between all molecules in the sample have a finite cutt-off 
radius of 15 A. In 3D simulatons, the contributions of sites outside this sphere are taken 
into account assuming an uniform distribution of sites and performing a simple integration. 
In our cuasi 2D system, the volume to integrate is that outside the cut-off sphere and within 
the volume of the slab. Appendix A includes this integral. 

b) The electrostatic interactions are calculated via the standard 3D Ewald sums with 
a large size box along the perpendicular to the bilayer slab and 2D periodic boundary 
conditions in the plane of the slab. The Ewald's sum term corresponding to the macroscopic 
electric field is discussed in the following section. 

IV. THE MACROSCOPIC ELECTRIC FIELD IN A QUASI - 2 D SAMPLE: 



Electrostatic forces have a far from negligible contribution to the self-assembly and final 
patterns found in soft matter systems. For bilayers the macroscopic electric field is given, in 
a first approximation, by the contribution of the surface charges of the macroscopic dielectric 



slab. There are very clear reviews where the first multipolar moment of these charges are 
taken into account, in 2D and 3D macroscopic samples [25, 27]. In particular, for an uniform 
dielectric slab oriented perpendicular to the z direction, the contribution of the uniformly 
distributed surface charges to the total energy of the macroscopic slab system, of volume V, 
is: 

T jmacrosc. 27r j,ir 2 

where Mz is the total dipole moment of the slab. Then the contribution of this term to the 
total force on every charge qi of the sample is: 

This approximation has been tested in monolayers [26, 27]. In the MD simulation of 
bilayers, and due to its geometry, the total dipole moment M^ is zero in a time average. 
Nevertheless, as in monolayers the forces p^<^crosc. ^^-^^ ^^^ ^^^ negligible by comparison 
with the rest of the interaction forces, and, as the electrostatic interactions play a key role 
in the dynamics and structural properties of the bilayers, a more accurate estimation of 
their macroscopic electric field is desirable. One possibility is to include higher multipolar 
moments of the surface charges, but the convergence of the electrostatic interactions using 
total multipole series is slow. Another possibility is to solve the Poisson-Boltzman equation 
with boundary conditions given by the surface charges [38], this method is extremely 
lenghtly for a numerical simulation, because the charged atoms are mobile and change their 
location in every time step. 

Nagle & Nagle [39, 40] recently reviewed the experimental data on the structure of lipid 
bilayers, in particular the distribution functions for the symmetrical components along 
the direction of the normal to the bilayer. They analize not only the peak positions of 
each distribution but their shape, pointing out that they are slightly asymmetrical (to the 
interior or exterior of the bilayer) and therefore a gaussian distribution model is just an 
useful first approximation when only the positions and widths are available. 

Here we propose a novel coarsed fit for the charge distribution of the different bilayer 
components (water and charged amphiphilics plus ions) using a superposition of gaussian 
distributions. We found that, in this way, the contribution of these charge distributions to 
the macroscopic electric field can be exactly calculated. The method is extremely simple to 
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implement in numerical simulations, and the spatial and temporal charge inhomogeneities 
are roughly taken into account. Our results can also be applied to biological membrane 
models. 

At each time step of the MD simulation we decompose our bilayer's charge distribution 
in four neutral slabs: two for the upper and lower water layers and other two for the neutral 
layers formed by the head plus ion charges. 

For each one of the four neutral slabs, instead of consider two planar surfaces with an 
uniform density of opposite charges, we propose two gaussian distribution along z, the 
perpendicular to the slabs, with the same opposite total charges and located the same 
relative distance (maintaining the slab width). The coarsed distribution of charges in the 
bilayer is then a linear superposition of gaussians: 

The macroscopic electric potential V{z) and the force field Ez{z) due to this type of 
charge distribution can be exactly solved. In Appendix B we include the analytically solved 
integrals (one of them is a new mathematical solution). The final result is: 

V{z) = -2v/2^E. -^ <1^ (exp(-^) - (^) Erf[i^]) . 

The obtained result is valid for any number of slabs, that as a function of time can 
change not only their position and width but also they can superpose. The practical 
implementation of this macroscopic field in our simulations (that is, how the ctj, Zi and qi 
values are taken at each time step) is explained in the following section. 
The contributions of this macroscopic field, as well as that of the external field Ueffih) 
in biological membranes, calculated for several bilayers samples are included in section 7. 
Units: density of charges [p] = -p; electrostatic potential [V] = 4- (for comparisson with 



experimental data 



,1, 



0.04803 ^'"'^V'^' = 14.399 [Volt] ); electric field [E] = -^ 



Here we have applied this exact calculation method of the macroscopic electric field to 
a symmetrical (along z) slab geometry, but it is also valid in an asymmetrical case, which 
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may imply a finite diff'erence of potential across the bilayer. 

The extension of this method ( a coarse grained representation of the macroscopic electric 
field via a superposition of gaussians) to other geometries is also strightforward, a spherical 
geometry, for example, would be useful for the study of miscelles. Its great advantage is 
that these representations can be analytically solved. Moreover, although the present work 
is not dedicated to the study of reaction fields, based on a coarse gaussian distribution of 
charges, can be also be easily applied to include a more realistic reaction field inside a cavity 
in non-uniform dielectrics. 

V. IMPLIMENTATION OF THE NUMERICAL SIMULATIONS: 

In our MD simulations of the bilayers of Fig. 1, the integration algorithms, time step 
and cut-off radius are essentially identical to those used on our bulk samples [33, 35], 
except for the periodic boundary conditions, that now are applied only in the xy plane 
of the bilayers, and that the calculations include now our proposed macroscopic electric 
field. The equations of motion of the rigid water molecules are integrated using the velocity 
Verlet algorithm for the atomic displacements and the Shake and Rattle algorithms for the 
constant bond length constraints on each molecule. The time step is of 1 fs., the sample is 
thermalized for 20 ps. and measured in the followings 100 ps. 

To maintain a constant temperature in our simulations of these bilayers, the Nose-Hoover 
chains method [41] (that we previously used for bulk samples) was disregarded because 
of the lenghly termalization of all the components, in many cases we obtained a different 
equilibrium temperature for each different molecular species. Instead, we chose to perform 
our MD simulations of the bilayers using the Berendsen algorithm [42], applying the 
equipartition theorem to each type of molecule and a strong coupling constant r = 0.5 At 
linking the average kinetic energy of each kind of molecules (amphiphilics, ions and water) 
to the desired kinetic energy of f/ceTo, with Tq = 300K. The Berendsen algorithm attains 
equipartition and equilibrium temperatures faster than the Nose-Hoover chains method, 
when applied to our mix of fiexible and rigid molecules. 

To include the macroscopic electric field term we need to determine the values of the ai, 
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Zi and qi parameters at each MD time step. For each one of the two slabs that simulate 
the charge distribution of chains' heads plus ions, we fit the Zi parameters of two gaussians, 
so as to reproduce the dipolar moment of the slab. Their Oi values are obtained from the 
corresponding charge distributions, with g^ = 1 for ions and gj = — 1 for chains' heads. 
Typical values of these variables, as well as the contribution of the external potential and 
the macroscopic electric field to the total forces on all molecules, are reported for each 
calculated case. 

VI. RESULTS: 



Here we present the results obtained for one sample of our simple NB film model. Our 
simulations were mainly performed to analyze the contribution of the macroscopic electric 
fields to the equilibrium structure and molecular dynamics of the bilayer. Elsewhere we 
will test the versatility of our approach by studying biological membrane models with and 
without ions solved in water, with and without the water layer, as a function of chain length, 
etc.. 





Figure 2: Final configuration of the HD sample with a=b=45A: (a) ac cross section, (b) ah cross 
section. 

In Ref. [14] a MD simulation with 64 all-atom sodium-dodecyl-sulfate SDS 
[Na~''CH3(CH2)ii(OS03)"] molecules and 1 to 6 water molecules per amphiphilic 
was presented. Here we use our simplified amphiphilic model. In section 4, Fig. 1 (h) 
showed the initial configuration of a sample consisting of 226 charged amphiphilics and 
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Figure 3: Atomic density profiles (number of atoms per A^) for the (a) HD and (6) LD head group 
samples. 

226 Na'^ ions solved in 365 water molecules, the bilayer is perpendicular to the 2 MD box 
axis, with a box size of a = 6 = 45. 5A, c =1000 A. The film is thermalized for 50 ps. and 
measured in the following 100 ps.. Fig. 2 shows the attained final configuration of this 
sample (HD high density), two other samples with a lower density of amphiphilic heads per 
unit area were calculated: a medium density MD (a = b = 46. 5A) and a LD low density 
sample (a = b = 48. OA). 

The width of our HD film is 38. 2A, of our MD film is 36. 8A and that of the LD one is 
35.1 A, as measured from the average distance between end tail CH^ groups. Fig. 3 shows 
the atomic density profiles in two cases. X-ray reflectivity measurements determined that 
common soap films (CBF) are several thousand angstroms wide, while a NBF is about 33 
A [15]. In ref. [43] a dual optical multiwave interferometer ( of A = 1.064/im) was used to 
determine the dynamics of gravity induced gradients in soap film thicknesses, which allows 
to detect variations of about lA in the thickness of the film. 

The three cases correspond to a glassy phase, as measured by their diffusion coefficients. 
For the HD sample the molecular centers of mass location show average oscillations with an 
amplitude of ~3 A and we measured an overall displacement of about 5 A in 100 ps., within 
the bilayer plane, therefore our measured diffusion coefficient is less than our measurement 
error (about lO^^cm^/sec). The experimental value of the diffusion coefficient for SDS, 
in a film of about 35A thickness, is 6 lO^^cm^/sec [7, 44] . In the MD and LD samples 
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Figure 4: Reorientational autocorrelation in the xy plane. 

we observe a collective rotation or large amplitud oscillation of the amphiphilic molecules 
around their corresponding center of mass and about an axis perpendicular to the bilayer 
plane. We measure a tilt angle of 26.2 deg. in our HD sample, 35.9 deg. in the MD and 
40.2 deg. in the LD sample. Fig. 4 shows the autocorrelation function of the xy component 
of the reorientational molecular motion, showing a disordered reorientation at HD and an 
oscillatory motion at LD. 

As in the previous bilayers models, the funcion p{z) is obtained from a coarse grained fit 
of the charge distribution of our MD sample using four charged slabs. The parameters that 

describe these slabs are for the LD film are: 

qi{e) Zi{A) cri{A) 

-1.0 3.485 3.344 
1.0 2.733 3.344 



slab 1 (head and ions) 
slab2 {water) 
slabs (water) 
slab 4 {head and ions) 



-0.241 
0.241 
0.241 

-0.241 

1.0 

-1.0 



1.196 0.433 

3.762 0.433 

-3.762 0.433 

-1.196 0.433 

-2.733 3.344 

-3.485 3.344 



and with them we calculate the functions included in Fig. 5 (a) . As this is a very thin NB 
film with about 1.6 water molecules per amphiphilic, water is highly packed and strongly 
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Figure 5: Profile of: (a) the macroscopic charge density p{z), electric field and potential (full line 
for total values, dotted line for water contribution), (b) our measured profile of macroscopic electric 
and total forces on molecules, as a function of the location of their centers of mass along z, in the 
HD, NB film. 

polarized due to the electric field generated by the charged head and ions. In Fig. 5 (a) 
we also included the strong contribution of the water layer to the calculated profiles of our 
samples. This behavior has a strong dependence on the content of water in the NBF and 
will be presented elsewhere. Our calculated total electrostatic potential is nearly constant 
at the core, with a negative value of about -1.69 e/A. Although not directly comparable, in 
Ref. [19] thicker NB films with up to 11.94 water molecules per amphiphilic were studied 
with an all-atom MD simulation and also a strong polarization of water was found at the 
interface with the head layers. Their measured electrostatic potential is also nearly constant 
at the core, but with a much smaller negative value of about -0.2 Volt= -0.0139e/A for the 
sample of 11.94 water molecules per amphiphilic. Unfortunately there is not experimental 
data available for comparison, but nevertheless, in Ref. [45] the origin of the short-range 
and strong repulsive force between two ionic surfactant layers is calculated as decaying 
exponentially with their distance. 

Fig. 5 (b) includes the contribution of macroscopic electric field E to our measured total 
forces on water, ions and chain molecules, as a function of the centers of mass location, 
averaged over a MD free trajectory of 100 ps. Units: [F] = kJ/mol/A. Due to the 
macroscopic electric field, the negatively charged amphiphilics tend to drift away of the 
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bilayer but this tendency is balanced by a contrary force on the movile ions, that remain 
within the water layer. When adding all molecule - molecule interactions, the total forces 
(an order of magnitude lower than those due to the macroscopic field) on the amphiphilics 
tend to maintain a bilayer structure with the water inside. 



VII. CONCLUSION: 

In this paper we studied a simple model of a amphiphilic bilayer, a Newton black film, 
and analyzed the key role of the electrostatic interactions in their self-assembly. The model 
is simple enough so a large number of molecules can be included in the MD samples, but 
also detailed enough so as to take into account molecular charge distributions, flexible 
amphiphilic molecules and a reliable model of water. All these properties are essential to 
obtain a reliable conclusion. 

The presented amphiphilic molecular model is extremely simple, but includes all 
main characteristics of these type of molecules: a strongly charged head atoms and a 
tail consisting in a semiflexible chain of hydrocarbon united atoms. As an example of 
its versatility we studied samples with a tail of 12 atoms and other of 18 atoms. This 
amphiphilic model also allows to study, in a simple way, the properties of bilayers formed 
by charged or neutral amphiphilics and with or without explicity including water molecules 
in the numerical simulations. 

As for the calculation of the electrostatic interactions, we also propose a novel and more 
accurate method to calculate the macroscopic electric held in cuasi 2D geometries, which 
can be easily included in any numerical calculation. The method, that essentially is a 
coarse grain flt of the macroscopic electric fleld beyond the dipole order approximation, 
was applied here to symmetrical bilayers (along the normal to the bilayer and with periodic 
boundary conditions in two dimensions), but their derivation is general and valid also for 
asymmetrical slab geometries. 

Lastly, we emphasize the relevance and utility of these simple bilayer models. They can 
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be applied to the systematic study of the physical properties of these bilayers, that strongly 
depend not only on 'external parameters', like surface tension and temperature, but also on 
the kind of guest molecules embedded in them (eg. proteins, etc.). In turn, the structure 
and dynamics of the embedded molecules strongly depend on their interactions with the 
bilayer and with the sorrounding water. 

These simple bilayers can also be extremely useful to model, for example, the synthesis 
of inorganic (ordered or disordered) materials via an organic agent, which are based on the 
use of these self-assembling molecules [9]. This is a recent and very fast growing research 
field of nanotechnologycal relevance, as are lithography, etching and molding devices at the 
nanoscopic scale. Another fast developing field is that of electronic sensors and nanodevices 
supported on lipid bilayers [4, 46]. 

Furthermore, the amphiphilic bilayers are hosts for the diffusion and/or nucleation of 
guest molecules of relevance in biological and/or ambient problems, these studies are cur- 
rently being performed. 

Lastly, as our model retains the fiexibility of the original amphiphilics, and the electro- 
static interactions are included, the approach is really useful to obtain 'realistic' solutions 
to the above mentioned problems as well as those related to electric fields and electrostatic 
properties. 
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A: Correction terms for the atom-atom interactions within the slab and beyond the 
cut-off radius 

In 3D samples, correction terms to the energy U^^{Rcut) and virial Vir'^'^{Rcut), due to 
the finite cut-off radius of the LJ site-site interactions, are calculated asuming an uniform 
distribution of sites beyond the cut-off radius Rcut- These terms are: 



U^'^'iRcut) = J^ sin edej^ ^ UlAt) 27rr2 dr = fea 



3 V -R, 



Re 



and 
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Vzr^'^iR.ut) = fo sin 9 dO /;j^^ (-r^^^) 27cr^ dr = Syrea^ 



Re 



Rcut 



In our 2D samples, instead, the integrals are over the volume outside the Rcut radius and 
inside the slab. In all of our cases the width of the slab Zgiab is larger than the diameter of 
the cut-off sphere {zsiah > '^Rcut), and the correction terms are then simply calculated. 

f/2^(i?e„*) = j; sin 9 d9 Jl^2f'^ UUr) 2vrr2 dr, 
where 9 is the azimutal angle measured from the perpendicular to the bilayer and Rmax (9) = 



2cos(0)- 

The final result is: 

U'^iRcut) = U'^'iRcut) - f eo-3 



30 I Zslab 



1 / 2cr 

4 I Zslab 



In a similar way we obtain: 
V^r^^{Rcut) = f;sm9d9j;^-^'\-r'-^)2nr^dr 



Vtr^''{Rcut)-^7rea' 



and for the component 
bilayer) : 



> 9 / \ 3 

2a \ 1 I 2a 

2 \ Zslab 



^ ' JU \^slab / ^ \Zslab / 

along z (used to control the pressure along the perpendicular to the 



Vir-riRcut) =j;sin9d9j;^- 
= ^Vtr^^'iRcut) 



W(_^at^^27rr2rfr 



fea^ 



2a 



2a 

^slab 



On the other hand, if Zsiab < 2Rcut , the calculated terms are: 

U'''{Rcut) = U^''{R 



^cut) OR ^ g ^^ ^slab 



2Rc. 



30 V R, 



Vir'^'^^Rcut) 
Virl^{Rcut) 






Vir^^iR, 



[jf^filVtr^^ (Rcut)- U'"" {Rcut 



4_ I 2a 

30 I Zslab 



4 yRcut 

9 



1 



2a 

^slab 



and 



B: The electrostatic macroscopic field model 

We propose a superposition of slabs, upper and lower surfaces with opposite charges, but 
with a gaussian distribution, of width a, along the depth of the slab. The charge distribution, 
electrostatic potential and electric field expresions are given in the following subsections, in 
order of increasing complexity. Case C.2 is the one we use for the reaction field in our MD 
simulations. 

The results showed here correspond to a symmetrical distribution of charges, in the 
bilayer, about the origin (i. e. identical charges, widths and location of the maximums 
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along ±z axis. But for an asymmetrical case the calculation is straightforward, and the 
obtained result should be very useful for studying these bilayers under an applied external 
electric field, for example. 

Bl: The macroscopic electric field of a single slab: 

First of all we review the well known case of an uniformly charged surface (p per unit 
area) located at z = 0. The potential field at z, due to this infinite, uniformly charged (xy) 
plane is: 

and the electric field is: 

m = -^= 2vrpft. 

For a homogeneous dielectric single slab (of width 2zi), the macroscopic electric field 
is calculated assuming an uniformly charged plane at Zi {qi per unit area) and another 
with opposite charge {—qi per unit area) located at -Zi. The potential field of this charge 
distribution is: 

V{z) = 2nqi {[V^Tz-z^r " l^Tz+^^r^ = -2nqi{\z - Zi\ - \z + Zi\), that is, 

-AirqiZi if z < -Zi 

^(^) — \ 4:7iqiZ if —Zi ^ z ^ Zi ■ 

AnqiZi if Zi> z 

v 
The electric field of this charge distribution is: 

E{z)=-4:7iqi, ii—Zi < z < Zi, and E{z) = if z < —Zi or z > Zi. 

Usually this result is given in terms of the the total dipole moment of the slab per unit 
area, which is M^ = 2qiZi, the macroscopic electric field then is E^ = — — M^ and the 
contribution of this charge distribution to the total energy of the system is ^M|. 

If instead of considering two uniformly opposite charged planes, we consider that the 
charges show a certain spread along z, we can analyze a model of two gaussian charge 
distributions along z, of total charge ig, (per unit area) each one, standard deviation ai and 
maximum located at ±2;irespectively, and the charge distribution across the slab becomes 
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now: 



pW = ^(-P(-^) --P(-^)) 



The total dipole moment of this system is again M^ = 2qiZi, but the potential V{z) will 
differ. For the sake of simplicity we take as equal the dispersion of both distribution and 
we obtain: 



V{z) = -2n JZ^ p{t) \z -t\dt = -2vr(/:^ p{t){z - t)dt + j^ p{t){t - z)dt). 

Replacing the p(t)function we first obtain: 

V{z) = -v^f ij:^ exp(-(^)(^ - t)dt - j:^ exp(-(^)(. - t)dt+ 

r exp(-(^)(t - z)dt - r exp(-(^)(t - z)dt ) , 
and finally the potential is: 

V{z)= -2v^ag,(exp(-(^)-exp(-(^)) 

-{z - z,) 27rg, {Erf[^] - Erf[^]) . 
From this function, the electric field across a single slab, centered at 2; = 0, is: 

When (7^0, the functions tend to those found in electrostatics textbooks. 



B2: Several neutral slabs model for the quasi-2D macroscopic electric field: 

The typical charge distribution of our samples of amphiphilic bilayers can be decomposed 
in a linear superposition of neutral slabs. Replacing for each slab the uniformly charged 
faces by gaussian distributions of charges along z, of total charge ±gj each one, width ctj 
and located at ±Zi, the preceding formulae in C.l can be easily extended to any number of 
slabs and the final results are: 

pW = E.7fcexp(-i^)) , 

V{z) = -2v/2^E. ^.g.(exp(-^) - (^) Erf[i^]) . 

The accumulated energy of this superposition of slabs is: 

w = UZ p{^)v{z)dz = - E. E, ?.?.^(/i + ^h). 

The integral Ji can be solved using the well known property of gaussians: the product of 
two gaussians is a third gaussian, in this case: 
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h = IZ exp(-(^)exp(-i^)rf. = ^exp(-,gg^); with 7 



V7 



J V ' ^v^^z T""j ;' 'j% <jj 



The second integral is a new mathematical solution and was solved with the aid of refs. 
[47, 48]. From the first one [47] we obtain: 

ha = jToo exp(— a(t — cf)erf{t)dt = y^er/( 'i^ )-, and from the second reference: 

hb = JZo t exp{-~a(t-c)'^)erf(t)dt = c/2a + ^^7f^ exp(-ac2 + Y^). Finally, and taking 

a = {aj/aiY and c = \)^J , we calculate: 

This formulae is extremely simple to include, and fast to calculate in a MD or MC 
simulation. The application to other geometries is also straightforward. These relationships 
are also valid to apply when working with systems where there is a finite difference of voltage 
through the slab, keeping in mind that the total charge of the slab is zero. 
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